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1. INTRODUCTION 

Numerical relativity represents the only currently viable method for obtaining 
solutions to Einstein's equations for highly dynamical and strong field sources 
of gravitational radiation. The most astrophysically interesting example is 
probably the final stages of binary inspiral and coalescence. Partly motivated 
by the prospect of observations with the next generation of gravitational wave 
detectors, a multi-institutional "Grand Challenge" effort is underway in the 
US aimed at solving the full Einstein equations numerically for coalescing black 
hole binaries. In addition to the tremendous computational resource demands 
of this problem, it has been realized for some time that, unfortunately, the 
standard 3+1 formulation of Einstein's theory as a Cauchy problem (cf. ||]) is 
somewhat deficient because of the difficulty of imposing boundary conditions 
which maintain numerical accuracy (and in some cases the physical correctness) 
of the solution. 

The problem of boundary conditions is most dramatically evident in the 
study of black hole spacetimes. Inside black holes spacelike slices either a) run 
into singularities causing termination of simulations, b) freeze their evolution 
necessitating the commitment of more and more computational resources to the 
astrophysically irrelevant black hole throat as the simulation progresses. An 
obvious solution to this problem is to excise the interior of the black hole from 
the computational domain. Since it is impossible to identify the event-horizon 
dynamically during the course of a simulation, a possible alternative is to use 
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the apparent horizon (which can be located on a single timeslice) and always 
lies inside the true horizon (assuming cosmic censorship holds). 

Numerical techniques based on the notion of causal differencing (cf. Sei- 
del/Suen in this volume and |3j) have been proposed for dealing with apparent 
horizon boundary conditions. In practice, it seems clear that the success of 
these techniques is crucially dependent on the form of the Einstein equations 
used{| ||. For spacetimes including gravitational radiation a purely hyperbolic 
evolution system is imperative because boundary conditions for the full set of 
constraint equations are not available on the apparent horizon. Furthermore, 
a purely hyperbolic evolution scheme with "simple" characteristics, where the 
only nonzero propagation speed is the speed of light, enables one to ignore 
entirely the spacetime inside the apparent horizon without concern that any 
unphysical (gauge) field quantities should be escaping the horizon. 

Outer boundary conditions for simulations on spacelike slices of asymptoti- 
cally flat spacetimes are another important issue for the computation of grav- 
itational waveforms. Since it is not feasible to simulate out to spatial infinity 
where there is no radiation, it is important to have boundary conditions that 
allow radiation to pass cleanly off the mesh. If an outgoing boundary con- 
dition is applied to the wrong variable, spurious radiation is produced which 
can contaminate the computed gravitational waveform. Additionally, for some 
problems it is necessary to put the outer boundary at such a small radius from 
the isolated source that backscatter of radiation off curvature is significant. 
This source of incoming radiation then needs to be built into the outer bound- 
ary conditions. The usual approach is to match the interior numerical solution 
onto perturbation theory for the exterior region §, |, |, §, @. Work is also 
underway to allow the connection of interior Cauchy solutions to exterior nu- 
merical solutions on characteristic hypersurfaces Jll]| . Both approaches benefit 
greatly from the use of a hyperbolic evolution scheme with simple characteris- 
tics for the interior solution. Outer boundary values can be assigned without 
the necessity of forming complicated admixtures of gauge and physical data. 

In this paper we motivate and discuss a remarkable new hyperbolic formula- 
tion of general relativity (lj, [l3) which may be thought of as a natural extension 
of the usual 3+1 split of spacetime. This formulation preserves complete spatial 
covariance by means of an arbitrary shift vector. The standard 3+1 treatment 
1 14, ||, is gauge covariant in this sense but not hyperbolic. Naturally, our 
formulation does require a condition on the time slicing to deal with the time- 
reparametrization invariance of the theory. This is physically intuitive; for 
example, we believe that a complete understanding of the generality of slicing 
conditions is a necessary first step towards addressing the problem of time in 
quantum gravity. We expect many other applications of this formulation. Early 
indications are that it will provide a powerful new approach to perturbation 
theory and approximation schemes for general relativity. 

The plan of this paper is as follows. First we will motivate the derivation of 
wave equations for general relativity by considering the vastly simpler case of 
a scalar wave and show how causal boundary conditions can be implemented. 
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We then turn our attention to electromagnetism to demonstrate the general 
procedure for gauge theories for producing wave equations. We then take gen- 
eral relativity in 3+1 form and apply the same method to obtain an explicitly 
hyperbolic formulation. This formulation is then written in first-order sym- 
metric hyperbolic form ideal for numerical implementation. Finally we discuss 
perturbative reductions of these equations and their use in outer boundary 
conditions and radiation extraction. 



2. BOUNDARY CONDITIONS FOR THE SCALAR WAVE EQUA- 
TION 

Consider the simple scalar wave equation in fiat space: 

□^ = 0. (1) 

Boundary conditions for this equation are straightforward to state and im- 
plement because the equation is manifestly hyperbolic. Since the equation is 
linear, it is clear how to impose outgoing wave boundary conditions at the edge 
of a numerical domain. It is also possible to employ inner "no boundary" or 
causal boundary conditions at the edge of an expanding null-surface (analo- 
gous to a black hole). Not surprisingly, since the equation is purely hyperbolic, 
stable and accurate solutions can be obtained by merely ignoring the causally 
disconnected region inside the null boundary. 

We have performed tests of this notion using the simple scalar wave equation 
(and straightforward nonlinear extensions) in 2+1 dimensions and Cartesian 
coordinates: 

d 2 tp d 2 ij) d 2 ip 

W = !h? + 'df ^' 
An arbitrary point in the grid is used as the origin of an expanding spherical 
wavefront that itself is used as the inner boundary. Outgoing wave boundary 
conditions are imposed at the rectangular outer boundary. Equation (||) is 
integrated using a predictor-corrector scheme. The second spatial derivatives 
are computed using centered differences except at the boundaries. At every 
timestep, the expanding inner boundary is located, boundary points identified 
and special finite difference operators constructed as shown in Figure [l| Enough 
points are used so that the 2nd-order error term can be set equal to the 2nd- 
order error in the centered difference scheme used in the mesh interior. 

We also simulate the effect of a shift-vector or grid velocity by allowing 
the coordinate system to be time dependent. This is accomplished by a high- 
order interpolation step following each time-step update. Since this shift-vector 
can be larger than the speed of propagation, it is possible for grid points to 
"emerge" from inside the horizon. These points are filled with data from out- 
side using high-order extrapolation. This code has been extensively tested 
with grid velocities up to 5 times the propagation speed. It is stable, second- 
order convergent, and equal in accuracy to a comparison code using standard 
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Fig. 1. — Schematic diagram of numerical differencing around a causal boundary. 



boundary conditions. Nonlinear source terms seem to present no problems for 
this scheme. A similar, equally successful, algorithm has been developed in 
the context of a flux-conservative first-order version of (g) and a Lax-Wendroff 
evolution scheme [|l5|. 

3. DEVELOPMENT OF HYPERBOLIC SYSTEM 

In the previous section we have demonstrated that hyperbolic wave equations 
are very amenable to imposition of causal boundary conditions. Here we dis- 
cuss the construction of analogous wave equations for considerably more com- 
plicated (and gauge dependent) theories. 
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3.1. Electromagnetism 

As a first example of a gauge theory, consider Maxwell's equations in flat space, 
written here in 3+1 form |]l6| . The dynamical equations are 

d t Ai = -Ek - (3) 

dtEi = -V j VjAi + ViV j Aj - 4vrJ, (4) 

where Ai is the vector potential, Ei is the electric field, Ji is the current source, 
V is a 3D flat space covariant derivative, and <j> is the gauge variable (scalar 
potential). These equations are supplemented by the initial value constraint 

V ] Ej = Airp (5) 

with p the charge density, and a gauge condition on Ai which requires the 
computation of <f>. To produce a wave equation, one approach is to take a 
time derivative of the Ai evolution equation and substitute the Ei evolution 
equation. To produce a D'Alembertian operator, it is necessary to apply, for 
example, a transversality condition on Ai which in turn imposes a radiation 
gauge condition on 0: VMj = — > V 4 Vi0 = — Anp (using the continuity 
equation). We have obtained a wave equation for Ai in the "Coulomb gauge." 
Alternatively, we could employ the Lorentz gauge: dt<f> + V 1 ^ = V p A^ = to 
obtain a wave equation for A^ = (— (f>,Ai) 

An alternative, gauge-covariant, approach is to take a time-derivative of the 
evolution equation for the electric field: 

d?E t = V,V 3 (-£, - V,0) - V^ji-Ei - Vi4>) - d t Ji (6) 

and use the constraint (|5|) to eliminate the first term yielding the wave equation 

UEi = 4irVip-d t Ji. (7) 

Interestingly, Ai doesn't appear in (Q) ; the dynamics of electromagnetism have 
been cleanly separated from the gauge-dependent evolution of the vector and 
scalar potentials. 

3.2. General Relativity 

Consider a globally hyperbolic manifold of topology Exi? with metric g^. A 
foliation of this spacetime is defined by a closed 1-form uj = V a t where t is 
this coordinate time function and io has normalization ||w|| = — N 2 with TV the 
lapse function. The four-dimensional line-element associated with may be 
decomposed in the general ADM|l4|] form as 

ds 2 = -N 2 dt 2 + gij (dx l + (3 l dt) {dx j + (3°dt), (8) 

where f3 l is the shift vector. It is convenient to introduce the non-coordinate 
co- frame, 

0° = dt, 9 l = dx l + ftdt (9) 
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with corresponding dual (convective) derivatives 

d = d/dt - i d/dx i , d, = d/dx\ (10) 
In this non-coordinate basis the ADM metric takes the simple form: 

ds 2 = -N 2 (e ) 2 + gij e i e j . (ii) 

Note that [<9o,d,] = {di(3 k )dk — Coi k dk, where the C's are the structure 
functions of the co-frame, d0 a = -\C^ a 9 p A 0f. 

Instead of using the usual time-congruence d/dt = do + f3 k dk which follows 
the spatial coordinates to define our evolution direction, wc define a more 
natural time derivative for evolution Q 

d o =d o +/3 k d k -C =d/dt-C p , (12) 

where Cp is the Lie derivative in a time slice £ along the shift vector. In 
combination with the lapse as N~ 1 do, this is the derivative with respect to 
proper time along the normal to E, and it always lies inside the light cone, 
in contrast to d/dt. In addition, it has the useful property that it commutes 
with the spatial coordinate derivatives, [do, di] = 0. This time-derivative is 
particularly appropriate to our form of the equations as we have subtracted out 
the momentum constraints which, in the Hamiltonian formulation, turn out to 
be generated by the shift-vector. The dynamical variables in the standard 3+1 
decomposition are the 3-metric gij and the extrinsic curvature of the slice £ as 
defined by the relation 

d o9l3 = —2NKij, (13) 
In four dimensions, we can write Einstein's equation as 

= k(T„„ - \g^T\). (14) 

Using the moving basis defined above, ( [l4|) can be split up into constraints and 
evolution equations. The time-time part of the Einstein equation leads to the 
Hamiltonian constraint 

Rl-R\ =-R-H 2 + K ij K ij (15) 

where H = K\ and R is the trace of the 3-dimensional Ricci tensor Rij (barred 
quantities are always spatial in our notation with indices running from 1-3). 
The time-space parts of Einstein's equation yield the momentum constraints: 

iV- 1 i? 0i = VjKj - V t H. (16) 

The purely spatial parts of Einstein's equation give us the evolution of the 
extrinsic curvature: 



Rij = - jjdoKij + HKij - 2KuK\ -^%VjN + Rij . ! 1 7 ) 
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The standard 3+1 formulation of general relativity consists of the evolution 
equations (O) and ( |l7| ) with initial data (gij,Kij) satisfying the constraints 



(15) and (jig). These equations are supplemented by equations for the sources 
(if any) and for the kinematical variables (3 l and N. The "slicing" equation for 
N often determined by a condition on H using the trace of (|l7|). 

To derive a wave equation for general relativity, one could, for example, follow 
the classic method of eliminating from the unwanted second derivatives of 
9p,v by using the full spacetime harmonic condition = jl7| [ll| ^9|. One 
would then obtain a non-geometric D'Alembertian g a ^d a dpg^ v . This procedure 
is analogous to using the Lorentz gauge in Maxwell's theory. 

Instead, we shall follow the spatially gauge-covariant analog of the procedure 
that led to (g) and (0) for Maxwell's equations. The spatial metric gij is 
analogous to Ai, the shift /3 to </>, and the extrinsic curvature Kij of £ to Ei. 
The lapse, on the other hand, is a quantity found only in time-reparametrization 
invariant theories. We take a time derivative of the equations of motion and 
subtract spatial gradients of the momentum constraints, thus obtaining a new 
quantity Qij 

doRij — ViRoj — Vji?io = Qij- (18) 
In terms of the dynamical gravitational variables, fiy may be expressed as 

fiy = NOKij + Jij + Sij , (19) 

where □ = —N~ 1 doN~ 1 do + V fc Vfe is the physical wave operator for arbitrary 
(3 k . It is constructed from second proper time-derivatives and second covariant 
spatial-derivatives. The source term is given by 

Jij = d (HK t] -2K z k K jk ) + (N~ 2 doN + H)V z y ] N 

-2N' 1 (\7 k N)V (l (NK k j) ) + 3(V k N)\J k K tJ (20) 
+N- 1 K, ij V k (NV k N) - 2\? {l (K j) k V k N) + N- 1 HV l V ] N 2 
+2N- 1 {V (l H){V, j) N 2 ) - 2NK k (l R ])k - 2NR kijm K km . 

and contains no second-derivatives of the extrinsic curvature. 
The slicing term is given by 

= -N-^VjidoN + N 2 H). (21) 

For fly to produce a wave equation, Sij must be equal to a functional involving 
fewer than second derivatives of Kij. The most obvious way to guarantee this 
is to use the harmonic condition (cf. p(| when (3 k = 0) 

d N + N 2 H = 0. (22) 

(This can easily be generalized by adding an ordinary well behaved function 
f{t, x) to the right hand side. The slicing generality compatible with hyper- 
bolic evolution schemes is the subject of Ref. Imposing ( p2|) for all time 
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Table I. — Possible evolution systems 




System (type) 


Equations 


Initial Data 


System I 


<)(,<!, . = —2NKij 


9ij 


(hyperbolic) 


NDKij = fly - Jij + VjVj-/(t, x) 


Kij, d Kij 




d N + N 2 H = f(t,x) 


N 


System II 


®o9ij = -2NKij 


9i i 


(Mixed hyperbolic/ 


N^Kij — Jij Sij 




elliptic) 


d h = -V k V k N + N(R + H 2 - gVRij) 


H = h(t,x) 



amounts to imposing an equation of motion for N. The shift is completely 
arbitrary: it can be given as a function of space and time or solved for on each 
timestep based on some condition on the evolved variables. This "System I" 
is purely hyperbolic. Initial data consists of 3-metric and extrinsic curvature 
satisfying the constraints and proper time-derivative of the extrinsic curvature 
satisfying (|l7|). The initial value for the lapse must also be specified. 

An alternative to the harmonic condition is to specify the trace of the ex- 
trinsic curvature as a known function for all time H = h(x,t). This also 
eliminates the second derivatives of unknown functions in Sij and provides a 
time-dependent elliptic equation, 

d h = -W k W k N + N(R + H 2 - g ij R i3 ), (23) 

to solve for N on each timestep. The shift vector is still freely specifiable. 
This "System II" is mixed hyperbolic/elliptic. It is possible to prove Jl^| using 
the doubly contracted Bianchi identity that Systems I and II are completely 
equivalent to Einstein's theory. These systems are summarized in Table [|. 

It is easily seen that for first-order perturbations of static backgrounds, the 
evolution equation for the 3-metric ( |l3| ) and the wave equation for the extrinsic 
curvature ( |l9| ) become completely decoupled. We will explore this idea further 
in the section on perturbative reduction. This situation is analogous to the 
separation of the equation for Ai from the wave equation for Ei we saw in 
linear electromagnetism. It is also consistent with physical intuition about 
the separation of transverse wave motion from longitudinal fields in general 
relativity. Locally, the 3-metric provides a background on top of which the 
extrinsic curvature propagates. 

3.3. First-order hyperbolic form for System I 

In order to reduce System I to hyperbolic form it is necessary to define some 
new variables. (Here we restrict ourselves to the vacuum case and the simple 
harmonic slicing condition f(t,x) — 0.) We introduce a, = TV^ViiV — the 
acceleration of the local Eulerian observers (those at rest in the time slices) — 
its derivatives do; = N~ 1 dodi and a,ji = Vj-a^ = a^, as well as time and space 
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derivatives of the extrinsic curvature 

doKa = NLij (24) 

and Mkij = V k Kij. System I can now be cast in flux-conservative first-order 
symmetric hyperbolic formp2|, p"3[ . The 49 unknowns of the first-order system 
are g u, N , Lij, Mkij, a-i, &ji and aoi, and the equations are (|l3]), (|22]), 



(|17|), Q2J, and 

d L %3 -NV k M ki:j = J t3 , (25) 

d Q M kij -NV k L i:j = N[a k L lJ +2M k(l m K j)m (26) 
+2if m (iM,-) fc m - 2K m (iM m j) k 
+2K m(l {K m 3) a k + a f) K m k - a m K j)k )}, 

d a i = -N{Ha i + M ik k ), (27) 

doaji-NVjaoi = Na k [2M {ij) k - M k ij (28) 
+2a {l K ]) k - a k K i:j ] + Naja 0i , 

d a 0t - NV k a kl = N[-R k l a k + a l {H 2 -2K kl K kl (29) 
+2a k a k + 2a k k ) + 2a k a k t + HM lk k - 2K kl M tkl ], 

where Jy is computed using (po|). To reduce this system to completely first- 
order form, the 3-dimensional Riemann curvature appearing in Jij is expressed 
in terms of the 3-dimensional Ricci curvature using the identity 

Rmijk = %9m\jRk]i + ^9i[kRj]m + R9m[k9j]i- (30) 

The 3D Ricci tensor is then eliminated using ( |l7| ) rewritten in terms of first- 
order variables: 

/-',, = Rij + L i:j - UK,, + 2K 2k K k J + a, a, + a,,. (31) 



The 4D Ricci tensor is computed from sources using (|14|). To demonstrate 
strictly first-order form it is also necessary to make Christoffel symbols part of 
the system by introducing a background metric as in Ref . |l^] . We stress again 
that the shift vector is not one of the unknown fields; the form of the equations 
is completely independent of /3 fe . 

With this form of the equations it is possible to read off the characteristic 
speeds of the different fields and verify one's physical expectations about the 
propagating degrees of freedom. Since there are no spatial derivatives on the 
left-hand-side of their evolution equations, we see that gij, Kij, N, and en all 
propagate with zero speed with respect to the Eulerian observers: they glide 
up the the normal to the foliation driven by the dynamical sources. Only 
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the derivatives of the extrinsic curvature (Ly, M^-fe) and the derivatives of the 
acceleration (a G j, ay) propagate with the speed of light. These quantities all 
appear in the components of the spacetime Riemann tensor and thus represent 
tidal fields. 

Along with A. Anderson (UNC) we have performed numerical tests of this 
form of the equations on the simple dynamical problem of even and odd-parity 
cylindrical waves. A Lax-Wendroff scheme is easily coded for the 27 equations 
necessary to describe this system with complete spatial gauge freedom. We find 
that a stable and accurate evolution can be computed which is comparable with 
that obtained by solving the usual 3+1 equations in fully harmonic coordinates 
(see (2^] for this version of the equations) . 

4. PERTURB ATIVE REDUCTION AND OUTER BOUNDARY 
CONDITIONS 

Here we sketch the reduction of System I for perturbations of the Schwarzschild 
metric. Full details are given in Ref. pp| . This reduction has helped us eluci- 
date many aspects of the full theory and provides new insight into the nature 
of gauge-invariant perturbation theory. It also allows us to define a frame- 
work for both radiation extraction and outer boundary conditions based on 
Schwarzschild perturbation theory. Such a framework can be used in conjunc- 
tion with numerical simulations. 

4.1. First-order perturbation theory of Schwarzschild 

For first order perturbations of static Schwarzschild, we make the following 
decomposition: 



9ij = .'/-.- + h H ( 32 ) 

Kij = + Kij (33) 

N = N + a (34) 

i = + v\ (35) 



Tildes denote background values. The background metric and lapse take their 
standard static Schwarzschild values and the background extrinsic curvature 
and shift are zero. Unless otherwise noted, covariant derivatives are with re- 
spect to the background metric. 

The evolution of the 3-metric ( |l3| ) reduces to 

dtgij = (36) 
d t K 3 = -2N Kij +2V (i « i ). (37) 

Notice that this equation is entirely gauge dependent: the arbitrary choice 
of shift v l translates into arbitrary distortion of metric perturbations. The 
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harmonic condition for the lapse splits into: 

d t N = (38) 
d t a = v l diN~N 2 n. (39) 

The wave equation for the extrinsic curvature ( |l9| ) reduces to 

j^dtKij - N\7 k \7 k K tJ = -4V (iK^ VfciV + N~ 1 n l3 V k NV k N + 3\7 k NV k K i:j 

+K i3 V k V k N - 2n k (i V 3) V k N - 27V" 1 n k (i V 3) NV ' k N + 2 K V l V j N 
+Ad {l nd j) N + 2N- 1 nV l NV j N - 2NR k{l K k 3) - 2NR kl3m n km . (40) 

Here n = K l i and i?y and Rij k i are background, spatial Ricci and Riemann 
tensors respectively. This equation is entirely decoupled from the evolution 
of the 3-metric perturbation (|37|). It could be directly evolved in an exterior 
region as a perturbative version of the full evolution equations. (Note that 
for flat space, N = 1, (|o|) reduces to DKy = which has radiative solutions 
corresponding to first time-derivatives of the usual transverse-traceless metric 
perturbations for gravitational waves.) 

Alternatively, one can perform a decomposition of (||o|) in terms of tensor 
spherical harmonics and produce scalar wave equations for the different £, m 
mode combinations. Here, for simplicity, we restrict attention to the slicing 
independent odd-parity perturbations. For odd-parity perturbations, Kij is 
decomposed with two tensor spherical harmonics and two amplitude functions. 
The component K rr j, is expressed in terms of the amplitude function a x and 
angular functions as 

«V0 = a x (t, r) sin 0d g Y im . (41) 

Taking the r-(f> component of ( fiTf ) and utilizing the <f> component of the mo- 
mentum constraint, a 1-dimensional scalar wave equation purely in terms of 
the amplitude function a x is formed: 

[<9 t 2 - (1 - 2M/r) 2 dl - (2/r)(l - 2M/r)d r - 2M/r 3 + 3M 2 /r 4 

+ (1 - 2M/r){£(e+l)/r 2 - 6M/r 3 )] a x {t,r) = 0. (42) 

To form the standard Regge- Wheeler equation for odd-parity perturbations of 
Schwarzschild (cf. ||j|) one takes a time-derivative of this equation using 

d t n r(j> = -V r V A> + NR r4> (43) 

which is the perturbative reduction of ( |l7| ) for odd-parity perturbations. Here 
the covariant derivatives are with respect to the perturbed background and the 
3D Ricci tensor is computed from the perturbed metric. The variable dta x 
satisfies the usual Regge- Wheeker equation. We note that no work has been 
required to construct spatial gauge invariants. These come "for free" in our 
spatially covariant wave-equation. 
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The situation for even parity perturbations is somewhat more involved be- 
cause the lapse perturbations ( |39j ) couple into our wave equation via the trace 
of Kij and the harmonic slicing condition at the extrinsic curvature level. The 
same basic procedure holds, however. Using a tensor spherical-harmonic de- 
composition and the radial component of the momentum constraint, coupled 
ID wave equations are formed for projections of K rr and of k. Connection to the 
Zerilli equation can be made by taking a time-derivative of these equations. The 
usual gauge invariant perturbation equations for Schwarzschild spacetime|23| 
are seen, not surprisingly, to represent curvature evolution. 



4.2. Radiation extraction and outer boundary conditions 



Perturbation theory has proven to be a powerful tool for extracting physical 
information from numerically generated spacetimes (cf. [[| [7| ||, ^ ) . The basic 
scheme is to match the full nonlinear interior solution to perturbation theory 
along the timelike cylinder representing the boundary of the computational 
domain. This idea is illustrated in Figure 2. 

The main steps in this scheme are to 1) construct perturbatively gauge- 
invariant quantities from evolved code variables, 2) propagate these gauge- 
invariants to large radius to remove near-zone effects, 3) use information from 
step 2 to construct code variables at the edge of the mesh, thus providing outer 
boundary conditions. The construction of gauge-invariants makes it possible to 
use the same extraction procedure in conjunction with numerical simulations 
using different choices of spatial gauges. This follows closely the conceptual 
picture for the calculation of gravitational radiation from isolated sources laid 
out by Thorne |24| . The calculation of the strong field and dynamical source 
is performed by the numerical simulation. The overlap/matching region is in 
the nondynamical near zone region (within a typical wavelength of the source). 
The goal is to compute waveforms in the wave zone beyond which the geometric 
optics approximation can be used to propagate the waves. In the procedure 
shown in Fig. g, the waveform is read off the perturbative variables at the 
outer boundary of the exterior evolution. Effects of backscatter off background 
curvature between the outer boundary of the interior nonlinear solution and 
the outer boundary of the exterior perturbative solution have been taken into 
account in both the waveform and the boundary conditions imposed on the 
interior solution. 

As should be clear from the discussion in the previous section, the new hy- 
perbolic formulation elucidates the process of attaching the standard 3+1 vari- 
ables gij, etc. onto perturbation theory. This subject is explored fully both 



for weak- field and Schwarzschild perturbation theory in Ref.jlO . In the weak 



field case, the exterior perturbative evolution can be done analytically. For 
Schwarzschild perturbations this requires a straightforward numerical integra- 
tion using the same coordinate time steps as the interior evolution. Boundary 
data for the exterior equations is computed via multipolar projections of the 
components of and Ly. The Schwarzschild mass is found from the ADM 
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Fig. 2. — Schematic diagram of a simulation with a solution to the full Einstein 
equations in the interior matched onto a perturbation theory solution in the exterior 
(shaded region). During the course of the evolution, boundary data for the pertur- 
bative evolution is read off from the nonlinear solution at the inner boundary. Data 
from the perturbative solution is used in turn to construct outer boundary data for 
the nonlinear simulation. Approximate asymptotic waveforms are read off the per- 
turbative solution at large radius. 



surface integral performed near the edge of the interior mesh. 

To produce boundary data for the interior simulation, the components of 
Kij and are reformed from the perturbative variables using the momentum 
constraint equations. For System I, the lapse is determined by the harmonic 
slicing condition which ties it to the trace of the extrinsic curvature. So the 
lapse at the outer boundary is set directly from the exterior evolution. (If 
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an elliptic slicing condition is used, then the lapse at the outer boundary will 
be known given boundary values for the extrinsic curvature and an imposed 
condition on (p3|).) The boundary condition on the 3-metric is determined using 
([l3|), the known boundary values for the extrinsic curvature and the lapse, and 
the chosen boundary condition on the shift vector components. 

As mentioned in the previous section, an alternative procedure for the ex- 
terior evolution is to simply integrate Eq. fl4C| ) in the exterior on a 3D finite 
difference grid. This can be accomplished either using a Cauchy or characteris- 
tic formulation of the equation and a spherical polar topology numerical mesh 
for computational efficiency Since the coordinate singularity at r = will not 
be part of the evolution domain, the usual difficulties with numerical instabili- 
ties will be avoided. In addition, it will be sufficient to perform adaptive mesh 
refinement, if desired, in only the radial direction. Imposition of boundary val- 
ues is trivial for the extrinsic curvature and proceeds exactly as above for the 
other variables. 
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